home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Source Code / Libraries / DCLAP 4j / SeqPups / apps / clustalw.src / trees.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-12-17  |  31.7 KB  |  1,152 lines  |  [TEXT/R*ch]
  1. /* Phyle of filogenetic tree calculating functions for CLUSTAL W */
  2. /* DES was here  FEB. 1994 */
  3.  
  4. #include <stdio.h>
  5. #include <string.h>
  6. #include <stdlib.h>
  7. #include <math.h>
  8. #include "clustalw.h"
  9. #include "dayhoff.h"    /* set correction for amino acid distances >= 75% */
  10.  
  11. /*
  12.  *   Prototypes
  13.  */
  14.  
  15. extern void *ckalloc(size_t);
  16. extern void ckfree(void *);
  17. extern void error(char *,...);  /* error reporting */
  18. extern int getint(char *, int, int, int);
  19. extern void get_path(char *, char *);
  20. extern FILE * open_output_file(char *, char *, char *, char *);
  21. /*random number routines in random.c */
  22. extern unsigned long addrand(unsigned long);
  23. extern void addrandinit(unsigned long);
  24. extern void   getstr(char *,char *);
  25.  
  26. void init_trees(void);
  27. void guide_tree(void);
  28. void phylogenetic_tree(void);
  29. void bootstrap_tree(void);
  30. void compare_tree(char **, char **, int *, int);
  31. void tree_gap_delete(void); /*flag all positions in alignment that have a gap */
  32. int dna_distance_matrix(FILE *);
  33. int prot_distance_matrix(FILE *);
  34. void nj_tree(char **, FILE *);
  35. void print_tree(char **, FILE *, int *);
  36. void print_phylip_tree(char **, FILE *, Boolean);
  37. void root_tree(char **, int);
  38. void distance_matrix_output(FILE *);
  39. void two_way_split(char **, FILE *, int, int, Boolean);
  40. Boolean transition(int,int);
  41.  
  42. /*
  43.  *   Global variables
  44.  */
  45.  
  46.  
  47. extern double **tmat;     /* general nxn array of reals; allocated from main */
  48.                           /* this is used as a distance matrix */
  49. extern Boolean dnaflag;   /* TRUE for DNA seqs; FALSE for proteins */
  50. extern Boolean tossgaps;  /* Ignore places in align. where ANY seq. has a gap*/
  51. extern Boolean kimura;    /* Use correction for multiple substitutions */
  52. extern Boolean output_tree_clustal;   /* clustal text output for trees */
  53. extern Boolean output_tree_phylip;    /* phylip nested parentheses format */
  54. extern Boolean output_tree_distances; /* phylip distance matrix */
  55. extern Boolean empty;                 /* any sequences in memory? */
  56. extern Boolean usemenu;   /* interactive (TRUE) or command line (FALSE) */
  57. extern int nseqs;
  58. extern int max_aa;
  59. extern int first_seq, last_seq;
  60. extern int *seqlen_array; /* the lengths of the sequences */
  61. extern char **seq_array;   /* the sequences */
  62. extern char **names;       /* the seq. names */
  63. extern char seqname[];        /* name of input file */
  64. extern char phylip_phy_tree_name[FILENAMELEN+1];
  65. extern FILE     *phylip_phy_tree_file;
  66.  
  67. static double     *av;
  68. static double     *left_branch, *right_branch;
  69. static double     *save_left_branch, *save_right_branch;
  70. static int    *boot_totals;
  71. static int     *kill;
  72. static int     ran_factor;
  73. static int     root_sequence;
  74. static int     *boot_positions;
  75. static FILE     *clustal_phy_tree_file;
  76. static FILE     *distances_phy_tree_file;
  77. static char     clustal_phy_tree_name[FILENAMELEN+1];
  78. static char     distances_phy_tree_name[FILENAMELEN+1];
  79. static Boolean     verbose;
  80. static char     *tree_gaps;
  81.                      /* array of weights; 1 for use this posn.; 0 don't */
  82.  
  83. extern int boot_ntrials;        /* number of bootstrap trials */
  84. extern unsigned int boot_ran_seed;    /* random number generator seed */
  85.  
  86. void phylogenetic_tree(void)
  87. /* 
  88.    Calculate a tree using the distances in the nseqs*nseqs array tmat.
  89.    This is the routine for getting the REAL trees after alignment.
  90. */
  91. {    char path[FILENAMELEN+1];
  92.     int i, j;
  93.     int overspill = 0;
  94.     int total_dists;
  95.     static char **standard_tree;
  96.     char lin2[10];
  97.  
  98.     if(empty) {
  99.         error("You must load an alignment first");
  100.         return;
  101.     }
  102.  
  103.     get_path(seqname,path);
  104.     
  105. if(output_tree_clustal)
  106.     if((clustal_phy_tree_file = open_output_file(
  107.         "\nEnter name for CLUSTAL    tree output file  ",path,
  108.         clustal_phy_tree_name,"nj")) == NULL) return;
  109.  
  110. if(output_tree_phylip)
  111.     if((phylip_phy_tree_file = open_output_file(
  112.         "\nEnter name for PHYLIP     tree output file  ",path,
  113.         phylip_phy_tree_name,"ph")) == NULL) return;
  114.  
  115. if(output_tree_distances)
  116.     if((distances_phy_tree_file = open_output_file(
  117.         "\nEnter name for distance matrix output file  ",path,
  118.         distances_phy_tree_name,"dst")) == NULL) return;
  119.  
  120.  
  121.     boot_positions = (int *)ckalloc( (MAXLEN+1) * sizeof (int) );
  122.  
  123.     for(j=1; j<=seqlen_array[first_seq]; ++j) 
  124.         boot_positions[j] = j;        
  125.  
  126.     if(output_tree_clustal) {
  127.         verbose = TRUE;     /* Turn on file output */
  128.         if(dnaflag)
  129.             overspill = dna_distance_matrix(clustal_phy_tree_file);
  130.         else 
  131.             overspill = prot_distance_matrix(clustal_phy_tree_file);
  132.     }
  133.  
  134.     if(output_tree_phylip) {
  135.         verbose = FALSE;     /* Turn off file output */
  136.         if(dnaflag)
  137.             overspill = dna_distance_matrix(phylip_phy_tree_file);
  138.         else 
  139.             overspill = prot_distance_matrix(phylip_phy_tree_file);
  140.     }
  141.  
  142.     if(output_tree_distances) {
  143.         verbose = FALSE;     /* Turn off file output */
  144.         if(dnaflag)
  145.             overspill = dna_distance_matrix(distances_phy_tree_file);
  146.         else 
  147.             overspill = prot_distance_matrix(distances_phy_tree_file);
  148.               distance_matrix_output(distances_phy_tree_file);
  149.     }
  150.  
  151. /* check if any distances overflowed the distance corrections */
  152.     if ( overspill > 0 ) {
  153.         total_dists = (nseqs*(nseqs-1))/2;
  154.         fprintf(stdout,"\n");
  155.         fprintf(stdout,"\n WARNING: %ld of the distances out of a total of %ld",
  156.         (long)overspill,(long)total_dists);
  157.         fprintf(stdout,"\n were out of range for the distance correction.");
  158.         fprintf(stdout,"\n This may not be fatal but you have been warned!");
  159.         fprintf(stdout,"\n");
  160.         fprintf(stdout,"\n SUGGESTIONS: 1) turn off the correction");
  161.         fprintf(stdout,"\n           or 2) remove the most distant sequences");
  162.         fprintf(stdout,"\n           or 3) use the PHYLIP package.");
  163.         fprintf(stdout,"\n\n");
  164.         if (usemenu) 
  165.             getstr("Press [RETURN] to continue",lin2);
  166.     }
  167.  
  168.     if(output_tree_clustal) verbose = TRUE;     /* Turn on file output */
  169.  
  170.     standard_tree   = (char **) ckalloc( (nseqs+1) * sizeof (char *) );
  171.     for(i=0; i<nseqs+1; i++) 
  172.         standard_tree[i]  = (char *) ckalloc( (nseqs+1) * sizeof(char) );
  173.  
  174.     if(output_tree_clustal || output_tree_phylip) 
  175.         nj_tree(standard_tree,clustal_phy_tree_file);
  176.  
  177.     if(output_tree_phylip) 
  178.         print_phylip_tree(standard_tree,phylip_phy_tree_file,FALSE);
  179.  
  180. /*
  181.     print_tree(standard_tree,phy_tree_file);
  182. */
  183.     ckfree((void *)tree_gaps);
  184.     ckfree((void *)boot_positions);
  185.     ckfree((void *)left_branch);
  186.     ckfree((void *)right_branch);
  187.     ckfree((void *)kill);
  188.     ckfree((void *)av);
  189.     for (i=1;i<nseqs+1;i++)
  190.         ckfree((void *)standard_tree[i]);
  191.     ckfree((void *)standard_tree);
  192.  
  193. if(output_tree_clustal) {
  194.     fclose(clustal_phy_tree_file);    
  195.     fprintf(stdout,"\nPhylogenetic tree file created:   [%s]\n",clustal_phy_tree_name);
  196. }
  197.  
  198. if(output_tree_phylip) {
  199.     fclose(phylip_phy_tree_file);    
  200.     fprintf(stdout,"\nPhylogenetic tree file created:   [%s]\n",phylip_phy_tree_name);
  201. }
  202.  
  203. if(output_tree_distances) {
  204.     fclose(distances_phy_tree_file);    
  205.     fprintf(stdout,"\nDistance matrix  file  created:   [%s]\n",distances_phy_tree_name);
  206. }
  207.  
  208.  
  209. }
  210.  
  211.  
  212.  
  213.  
  214.  
  215. Boolean transition(int base1, int base2) /* TRUE if transition; else FALSE */
  216. /* 
  217.  
  218.    assumes that the bases of DNA sequences have been translated as
  219.    a,A = 0;   c,C = 1;   g,G = 2;   t,T,u,U = 3;  N = 4;  
  220.  
  221.    A <--> G  and  T <--> C  are transitions;  all others are transversions.
  222.  
  223. */
  224. {
  225.     if( ((base1 == 0) && (base2 == 2)) || ((base1 == 2) && (base2 == 0)) )
  226.         return TRUE;                                     /* A <--> G */
  227.     if( ((base1 == 3) && (base2 == 1)) || ((base1 == 1) && (base2 == 3)) )
  228.         return TRUE;                                     /* T <--> C */
  229.     return FALSE;
  230. }
  231.  
  232.  
  233. void tree_gap_delete(void)   /* flag all positions in alignment that have a gap */
  234. {              /* in ANY sequence */
  235.     int seqn, posn;
  236.  
  237.     tree_gaps = (char *)ckalloc( (MAXLEN+1) * sizeof (char) );
  238.         
  239.     for(posn=1; posn<=seqlen_array[first_seq]; ++posn) {
  240.         tree_gaps[posn] = 0;
  241.          for(seqn=1; seqn<=last_seq-first_seq+1; ++seqn)  {
  242.             if((seq_array[seqn+first_seq-1][posn] < 0) ||
  243.                (seq_array[seqn+first_seq-1][posn] > max_aa)) {
  244.                tree_gaps[posn] = 1;
  245.                 break;
  246.             }
  247.         }
  248.     }
  249. }
  250.  
  251. void distance_matrix_output(FILE *ofile)
  252. {
  253.     int i,j;
  254.     
  255.     fprintf(ofile,"%6d",(pint)last_seq-first_seq+1);
  256.     for(i=1;i<=last_seq-first_seq+1;i++) {
  257.         fprintf(ofile,"\n%-10s ",names[i]);
  258.         for(j=1;j<=last_seq-first_seq+1;j++) {
  259.             fprintf(ofile,"%6.3f ",tmat[i][j]);
  260.             if(j % 8 == 0) {
  261.                 fprintf(ofile,"\n"); 
  262.                 if(j != last_seq-first_seq+1 ) fprintf(ofile,"          ");
  263.             }
  264.         }
  265.     }
  266. }
  267.  
  268.  
  269.  
  270. void nj_tree(char **tree_description, FILE *tree)
  271. {
  272.     register int i;
  273.     int l[4],nude,k;
  274.     int nc,mini,minj,j,j1,ii,jj;
  275.     double fnseqs,fnseqs2,sumd;
  276.     double diq,djq,dij,d2r,dr,dio,djo,da;
  277.     double tmin,total,dmin;
  278.     double bi,bj,b1,b2,b3,branch[4];
  279.     int typei,typej,type[4];             /* 0 = node; 1 = OTU */
  280.     
  281.     fnseqs = (double)last_seq-first_seq+1;
  282.  
  283. /*********************** First initialisation ***************************/
  284.     
  285.     if(verbose)  {
  286.         fprintf(tree,"\n\n\t\t\tNeighbor-joining Method\n");
  287.         fprintf(tree,"\n Saitou, N. and Nei, M. (1987)");
  288.         fprintf(tree," The Neighbor-joining Method:");
  289.         fprintf(tree,"\n A New Method for Reconstructing Phylogenetic Trees.");
  290.         fprintf(tree,"\n Mol. Biol. Evol., 4(4), 406-425\n");
  291.         fprintf(tree,"\n\n This is an UNROOTED tree\n");
  292.         fprintf(tree,"\n Numbers in parentheses are branch lengths\n\n");
  293.     }    
  294.  
  295.     mini = minj = 0;
  296.  
  297.     left_branch     = (double *) ckalloc( (nseqs+2) * sizeof (double)   );
  298.     right_branch    = (double *) ckalloc( (nseqs+2) * sizeof (double)   );
  299.     kill         = (int *) ckalloc( (nseqs+1) * sizeof (int) );
  300.     av           = (double *) ckalloc( (nseqs+1) * sizeof (double)   );
  301.  
  302.     for(i=1;i<=last_seq-first_seq+1;++i) 
  303.         {
  304.         tmat[i][i] = av[i] = 0.0;
  305.         kill[i] = 0;
  306.         }
  307.  
  308. /*********************** Enter The Main Cycle ***************************/
  309.  
  310.  /*    for(nc=1; nc<=(last_seq-first_seq+1-3); ++nc) {  */                /**start main cycle**/
  311.     for(nc=1; nc<=(last_seq-first_seq+1-3); ++nc) {
  312.         sumd = 0.0;
  313.         for(j=2; j<=last_seq-first_seq+1; ++j)
  314.             for(i=1; i<j; ++i) {
  315.                 tmat[j][i] = tmat[i][j];
  316.                 sumd = sumd + tmat[i][j];
  317.             }
  318.  
  319.         tmin = 99999.0;
  320.  
  321. /*.................compute SMATij values and find the smallest one ........*/
  322.  
  323.         for(jj=2; jj<=last_seq-first_seq+1; ++jj) 
  324.             if(kill[jj] != 1) 
  325.                 for(ii=1; ii<jj; ++ii)
  326.                     if(kill[ii] != 1) {
  327.                         diq = djq = 0.0;
  328.  
  329.                         for(i=1; i<=last_seq-first_seq+1; ++i) {
  330.                             diq = diq + tmat[i][ii];
  331.                             djq = djq + tmat[i][jj];
  332.                         }
  333.  
  334.                         dij = tmat[ii][jj];
  335.                         d2r = diq + djq - (2.0*dij);
  336.                         dr  = sumd - dij -d2r;
  337.                         fnseqs2 = fnseqs - 2.0;
  338.                             total= d2r+ fnseqs2*dij +dr*2.0;
  339.                         total= total / (2.0*fnseqs2);
  340.  
  341.                         if(total < tmin) {
  342.                             tmin = total;
  343.                             mini = ii;
  344.                             minj = jj;
  345.                         }
  346.                     }
  347.         
  348.  
  349. /*.................compute branch lengths and print the results ........*/
  350.  
  351.  
  352.         dio = djo = 0.0;
  353.         for(i=1; i<=last_seq-first_seq+1; ++i) {
  354.             dio = dio + tmat[i][mini];
  355.             djo = djo + tmat[i][minj];
  356.         }
  357.  
  358.         dmin = tmat[mini][minj];
  359.         dio = (dio - dmin) / fnseqs2;
  360.         djo = (djo - dmin) / fnseqs2;
  361.         bi = (dmin + dio - djo) * 0.5;
  362.         bj = dmin - bi;
  363.         bi = bi - av[mini];
  364.         bj = bj - av[minj];
  365.  
  366.         if( av[mini] > 0.0 )
  367.             typei = 0;
  368.         else
  369.             typei = 1;
  370.         if( av[minj] > 0.0 )
  371.             typej = 0;
  372.         else
  373.             typej = 1;
  374.  
  375.         if(verbose) 
  376.              fprintf(tree,"\n Cycle%4d     = ",(pint)nc);
  377.  
  378. /* 
  379.    set negative branch lengths to zero.  Also set any tiny positive
  380.    branch lengths to zero.
  381. */        if( fabs(bi) < 0.0001) bi = 0.0;
  382.         if( fabs(bj) < 0.0001) bj = 0.0;
  383.  
  384.             if(verbose) {
  385.             if(typei == 0) 
  386.             fprintf(tree,"Node:%4d (%9.5f) joins ",(pint)mini,bi);
  387.             else 
  388.             fprintf(tree," SEQ:%4d (%9.5f) joins ",(pint)mini,bi);
  389.  
  390.             if(typej == 0) 
  391.             fprintf(tree,"Node:%4d (%9.5f)",(pint)minj,bj);
  392.             else 
  393.             fprintf(tree," SEQ:%4d (%9.5f)",(pint)minj,bj);
  394.  
  395.             fprintf(tree,"\n");
  396.             }    
  397.  
  398.  
  399.             left_branch[nc] = bi;
  400.             right_branch[nc] = bj;
  401.  
  402.         for(i=1; i<=last_seq-first_seq+1; i++)
  403.             tree_description[nc][i] = 0;
  404.  
  405.              if(typei == 0) { 
  406.             for(i=nc-1; i>=1; i--)
  407.                 if(tree_description[i][mini] == 1) {
  408.                     for(j=1; j<=last_seq-first_seq+1; j++)  
  409.                          if(tree_description[i][j] == 1)
  410.                             tree_description[nc][j] = 1;
  411.                     break;
  412.                 }
  413.         }
  414.         else
  415.             tree_description[nc][mini] = 1;
  416.  
  417.         if(typej == 0) {
  418.             for(i=nc-1; i>=1; i--) 
  419.                 if(tree_description[i][minj] == 1) {
  420.                     for(j=1; j<=last_seq-first_seq+1; j++)  
  421.                          if(tree_description[i][j] == 1)
  422.                             tree_description[nc][j] = 1;
  423.                     break;
  424.                 }
  425.         }
  426.         else
  427.             tree_description[nc][minj] = 1;
  428.             
  429.  
  430. /* 
  431.    Here is where the -0.00005 branch lengths come from for 3 or more
  432.    identical seqs.
  433. */
  434. /*        if(dmin <= 0.0) dmin = 0.0001; */
  435.                 if(dmin <= 0.0) dmin = 0.000001;
  436.         av[mini] = dmin * 0.5;
  437.  
  438. /*........................Re-initialisation................................*/
  439.  
  440.         fnseqs = fnseqs - 1.0;
  441.         kill[minj] = 1;
  442.  
  443.         for(j=1; j<=last_seq-first_seq+1; ++j) 
  444.             if( kill[j] != 1 ) {
  445.                 da = ( tmat[mini][j] + tmat[minj][j] ) * 0.5;
  446.                 if( (mini - j) < 0 ) 
  447.                     tmat[mini][j] = da;
  448.                 if( (mini - j) > 0)
  449.                     tmat[j][mini] = da;
  450.             }
  451.  
  452.         for(j=1; j<=last_seq-first_seq+1; ++j)
  453.             tmat[minj][j] = tmat[j][minj] = 0.0;
  454.  
  455.  
  456. /****/    }                        /**end main cycle**/
  457.  
  458. /******************************Last Cycle (3 Seqs. left)********************/
  459.  
  460.     nude = 1;
  461.  
  462.     for(i=1; i<=last_seq-first_seq+1; ++i)
  463.         if( kill[i] != 1 ) {
  464.             l[nude] = i;
  465.             nude = nude + 1;
  466.         }
  467.  
  468.     b1 = (tmat[l[1]][l[2]] + tmat[l[1]][l[3]] - tmat[l[2]][l[3]]) * 0.5;
  469.     b2 =  tmat[l[1]][l[2]] - b1;
  470.     b3 =  tmat[l[1]][l[3]] - b1;
  471.  
  472.     branch[1] = b1 - av[l[1]];
  473.     branch[2] = b2 - av[l[2]];
  474.     branch[3] = b3 - av[l[3]];
  475.  
  476. /* Reset tiny negative and positive branch lengths to zero */
  477.     if( fabs(branch[1]) < 0.0001) branch[1] = 0.0;
  478.     if( fabs(branch[2]) < 0.0001) branch[2] = 0.0;
  479.     if( fabs(branch[3]) < 0.0001) branch[3] = 0.0;
  480.  
  481.     left_branch[last_seq-first_seq+1-2] = branch[1];
  482.     left_branch[last_seq-first_seq+1-1] = branch[2];
  483.     left_branch[last_seq-first_seq+1]   = branch[3];
  484.  
  485.     for(i=1; i<=last_seq-first_seq+1; i++)
  486.         tree_description[last_seq-first_seq+1-2][i] = 0;
  487.  
  488.     if(verbose)
  489.         fprintf(tree,"\n Cycle%4d (Last cycle, trichotomy):\n",(pint)nc);
  490.  
  491.     for(i=1; i<=3; ++i) {
  492.        if( av[l[i]] > 0.0) {
  493.               if(verbose)
  494.                   fprintf(tree,"\n\t\t Node:%4d (%9.5f) ",(pint)l[i],branch[i]);
  495.         for(k=last_seq-first_seq+1-3; k>=1; k--)
  496.             if(tree_description[k][l[i]] == 1) {
  497.                 for(j=1; j<=last_seq-first_seq+1; j++)
  498.                      if(tree_description[k][j] == 1)
  499.                         tree_description[last_seq-first_seq+1-2][j] = i;
  500.                 break;
  501.             }
  502.        }
  503.        else  {
  504.               if(verbose)
  505.                fprintf(tree,"\n\t\t  SEQ:%4d (%9.5f) ",(pint)l[i],branch[i]);
  506.         tree_description[last_seq-first_seq+1-2][l[i]] = i;
  507.        }
  508.        if(i < 3) {
  509.               if(verbose)
  510.                 fprintf(tree,"joins");
  511.        }
  512.     }
  513.  
  514.     if(verbose)
  515.         fprintf(tree,"\n");
  516.  
  517. }
  518.  
  519.  
  520.  
  521.  
  522. void bootstrap_tree(void)
  523. {
  524.     int i,j,ranno;
  525.     int k,l,m,p;
  526.     unsigned int num;
  527.     char lin2[MAXLINE],path[MAXLINE+1];
  528.         char dummy[10];
  529.     static char **sample_tree;
  530.     static char **standard_tree;
  531.     int total_dists, overspill, total_overspill = 0;
  532.     int nfails = 0;
  533.  
  534.     if(empty) {
  535.         error("You must load an alignment first");
  536.         return;
  537.     }
  538.  
  539.     if(!output_tree_clustal && !output_tree_phylip) {
  540.         error("You must select either clustal or phylip tree output format");
  541.         return;
  542.     }
  543.     get_path(seqname, path);
  544.     
  545.     if (output_tree_clustal)
  546.         if((clustal_phy_tree_file = open_output_file(
  547.         "\nEnter name for bootstrap output file  ",path,
  548.         clustal_phy_tree_name,"njb")) == NULL) return;
  549.  
  550.     if (output_tree_phylip)
  551.         if((phylip_phy_tree_file = open_output_file(
  552.         "\nEnter name for bootstrap output file  ",path,
  553.         phylip_phy_tree_name,"phb")) == NULL) return;
  554.  
  555.     boot_totals    = (int *)ckalloc( (nseqs+1) * sizeof (int) );
  556.     for(i=0;i<nseqs+1;i++)
  557.         boot_totals[i]=0;
  558.         
  559.     boot_positions = (int *)ckalloc( (MAXLEN+1) * sizeof (int) );
  560.  
  561.     for(j=1; j<=seqlen_array[first_seq]; ++j)  /* First select all positions for */
  562.         boot_positions[j] = j;       /* the "standard" tree */
  563.  
  564.     if(output_tree_clustal) {
  565.         verbose = TRUE;     /* Turn on file output */
  566.         if(dnaflag)
  567.             overspill = dna_distance_matrix(clustal_phy_tree_file);
  568.         else 
  569.             overspill = prot_distance_matrix(clustal_phy_tree_file);
  570.     }
  571.  
  572.     if(output_tree_phylip) {
  573.         verbose = FALSE;     /* Turn off file output */
  574.         if(dnaflag)
  575.             overspill = dna_distance_matrix(phylip_phy_tree_file);
  576.         else 
  577.             overspill = prot_distance_matrix(phylip_phy_tree_file);
  578.     }
  579.  
  580. /* check if any distances overflowed the distance corrections */
  581.     if ( overspill > 0 ) {
  582.         total_dists = (nseqs*(nseqs-1))/2;
  583.         fprintf(stdout,"\n");
  584.         fprintf(stdout,"\n WARNING: %ld of the distances out of a total of %ld",
  585.         (pint)overspill,(pint)total_dists);
  586.         fprintf(stdout,"\n were out of range for the distance correction.");
  587.         fprintf(stdout,"\n This may not be fatal but you have been warned!");
  588.         fprintf(stdout,"\n");
  589.         fprintf(stdout,"\n SUGGESTIONS: 1) turn off the correction");
  590.         fprintf(stdout,"\n           or 2) remove the most distant sequences");
  591.         fprintf(stdout,"\n           or 3) use the PHYLIP package.");
  592.         fprintf(stdout,"\n\n");
  593.         if (usemenu) 
  594.             getstr("Press [RETURN] to continue",dummy);
  595.     }
  596.  
  597.     ckfree((void *)tree_gaps);
  598.  
  599.     if (output_tree_clustal) verbose = TRUE;   /* Turn on screen output */
  600.  
  601.     standard_tree   = (char **) ckalloc( (nseqs+1) * sizeof (char *) );
  602.     for(i=0; i<nseqs+1; i++) 
  603.         standard_tree[i]   = (char *) ckalloc( (nseqs+1) * sizeof(char) );
  604.  
  605. /* compute the standard tree */
  606.  
  607.     if(output_tree_clustal || output_tree_phylip)
  608.         nj_tree(standard_tree,clustal_phy_tree_file);
  609.  
  610.     if (output_tree_clustal)
  611.         fprintf(clustal_phy_tree_file,"\n\n\t\t\tBootstrap Confidence Limits\n\n");
  612.  
  613. /* save the left_branch and right_branch for phylip output */
  614.     save_left_branch = (double *) ckalloc( (nseqs+2) * sizeof (double)   );
  615.     save_right_branch = (double *) ckalloc( (nseqs+2) * sizeof (double)   );
  616.     for (i=1;i<=nseqs;i++) {
  617.         save_left_branch[i] = left_branch[i];
  618.         save_right_branch[i] = right_branch[i];
  619.     }
  620. /*  
  621.   The next line is a fossil from the days of using the cc ran()
  622.     ran_factor = RAND_MAX / seqlen_array[first_seq]; 
  623. */
  624.  
  625.     if(usemenu) 
  626.            boot_ran_seed = 
  627. getint("\n\nEnter seed no. for random number generator ",1,1000,boot_ran_seed);
  628.  
  629. /* do not use the native cc ran()
  630.     srand(boot_ran_seed);
  631. */
  632.            addrandinit((unsigned long) boot_ran_seed);
  633.  
  634.     if (output_tree_clustal)
  635.         fprintf(clustal_phy_tree_file,"\n Random number generator seed = %7u\n",
  636.         boot_ran_seed);
  637.  
  638.     if(usemenu) 
  639.           boot_ntrials = 
  640. getint("\n\nEnter number of bootstrap trials ",1,10000,boot_ntrials);
  641.  
  642.     if (output_tree_clustal) {
  643.           fprintf(clustal_phy_tree_file,"\n Number of bootstrap trials   = %7d\n",
  644.     (pint)boot_ntrials);
  645.  
  646.         fprintf(clustal_phy_tree_file,
  647.         "\n\n Diagrammatic representation of the above tree: \n");
  648.         fprintf(clustal_phy_tree_file,"\n Each row represents 1 tree cycle;");
  649.         fprintf(clustal_phy_tree_file," defining 2 groups.\n");
  650.         fprintf(clustal_phy_tree_file,"\n Each column is 1 sequence; ");
  651.         fprintf(clustal_phy_tree_file,"the stars in each line show 1 group; ");
  652.         fprintf(clustal_phy_tree_file,"\n the dots show the other\n");
  653.         fprintf(clustal_phy_tree_file,"\n Numbers show occurences in bootstrap samples.");
  654.     }
  655. /*
  656.     print_tree(standard_tree, clustal_phy_tree_file, boot_totals);
  657. */
  658.     verbose = FALSE;                   /* Turn OFF screen output */
  659.  
  660.     ckfree((void *)left_branch);
  661.     ckfree((void *)right_branch);
  662.     ckfree((void *)kill);
  663.     ckfree((void *)av);
  664.  
  665.     sample_tree   = (char **) ckalloc( (nseqs+1) * sizeof (char *) );
  666.     for(i=0; i<nseqs+1; i++) 
  667.         sample_tree[i]   = (char *) ckalloc( (nseqs+1) * sizeof(char) );
  668.  
  669.     fprintf(stdout,"\n\nEach dot represents 10 trials\n\n");
  670.         total_overspill = 0;
  671.     nfails = 0;
  672.     for(i=1; i<=boot_ntrials; ++i) {
  673.         for(j=1; j<=seqlen_array[first_seq]; ++j) { /* select alignment */
  674.                                 /* positions for */
  675.             ranno = addrand( (unsigned long) seqlen_array[1]) + 1;
  676.             boot_positions[j] = ranno;         /* bootstrap sample */
  677.         }
  678.         if(output_tree_clustal) {
  679.             if(dnaflag)
  680.                 overspill = dna_distance_matrix(clustal_phy_tree_file);
  681.             else 
  682.                 overspill = prot_distance_matrix(clustal_phy_tree_file);
  683.         }
  684.     
  685.         if(output_tree_phylip) {
  686.             if(dnaflag)
  687.                 overspill = dna_distance_matrix(phylip_phy_tree_file);
  688.             else 
  689.                 overspill = prot_distance_matrix(phylip_phy_tree_file);
  690.         }
  691.  
  692.         if( overspill > 0) {
  693.             total_overspill = total_overspill + overspill;
  694.             nfails++;
  695.         }            
  696.  
  697.         ckfree((void *)tree_gaps);
  698.  
  699.         if(output_tree_clustal || output_tree_phylip) 
  700.             nj_tree(sample_tree,clustal_phy_tree_file);
  701.  
  702.         compare_tree(standard_tree, sample_tree, boot_totals, last_seq-first_seq+1);
  703.         if(i % 10  == 0) fprintf(stdout,".");
  704.         if(i % 100 == 0) fprintf(stdout,"\n");
  705.     }
  706.  
  707. /* check if any distances overflowed the distance corrections */
  708.     if ( nfails > 0 ) {
  709.         total_dists = (nseqs*(nseqs-1))/2;
  710.         fprintf(stdout,"\n");
  711.         fprintf(stdout,"\n WARNING: %ld of the distances out of a total of %ld times %ld",
  712.         (long)total_overspill,(long)total_dists,(long)boot_ntrials);
  713.         fprintf(stdout,"\n were out of range for the distance correction.");
  714.         fprintf(stdout,"\n This affected %d out of %d bootstrap trials.",
  715.         (pint)nfails,(pint)boot_ntrials);
  716.         fprintf(stdout,"\n This may not be fatal but you have been warned!");
  717.         fprintf(stdout,"\n");
  718.         fprintf(stdout,"\n SUGGESTIONS: 1) turn off the correction");
  719.         fprintf(stdout,"\n           or 2) remove the most distant sequences");
  720.         fprintf(stdout,"\n           or 3) use the PHYLIP package.");
  721.         fprintf(stdout,"\n\n");
  722.         if (usemenu) 
  723.             getstr("Press [RETURN] to continue",dummy);
  724.     }
  725.  
  726.  
  727.     ckfree((void *)boot_positions);
  728.  
  729.     for (i=1;i<nseqs+1;i++)
  730.         ckfree((void *)sample_tree[i]);
  731.     ckfree((void *)sample_tree);
  732. /*
  733.     fprintf(clustal_phy_tree_file,"\n\n Bootstrap totals for each group\n");
  734. */
  735.     if (output_tree_clustal)
  736.         print_tree(standard_tree, clustal_phy_tree_file, boot_totals);
  737.  
  738.     if(output_tree_phylip) {
  739.         for (i=1;i<=nseqs;i++) {
  740.             left_branch[i] = save_left_branch[i];
  741.             right_branch[i] = save_right_branch[i];
  742.         }
  743.         print_phylip_tree(standard_tree,phylip_phy_tree_file,
  744.                                  TRUE);
  745.     }
  746.  
  747.     ckfree((void *)boot_totals);
  748.     ckfree((void *)save_left_branch);
  749.     ckfree((void *)save_right_branch);
  750.  
  751.     for (i=1;i<nseqs+1;i++)
  752.         ckfree((void *)standard_tree[i]);
  753.     ckfree((void *)standard_tree);
  754.  
  755.     if (output_tree_clustal)
  756.         fclose(clustal_phy_tree_file);
  757.  
  758.     if (output_tree_phylip)
  759.         fclose(phylip_phy_tree_file);
  760.  
  761.     if (output_tree_clustal)
  762.         fprintf(stdout,"\n\nBootstrap output file completed       [%s]"
  763.         ,clustal_phy_tree_name);
  764.     if (output_tree_phylip)
  765.         fprintf(stdout,"\n\nBootstrap output file completed       [%s]"
  766.         ,phylip_phy_tree_name);
  767. }
  768.  
  769.  
  770. void compare_tree(char **tree1, char **tree2, int *hits, int n)
  771. {    
  772.     int i,j,k,l;
  773.     int nhits1, nhits2;
  774.  
  775.     for(i=1; i<=n-3; i++)  {
  776.         for(j=1; j<=n-3; j++)  {
  777.             nhits1 = 0;
  778.             nhits2 = 0;
  779.             for(k=1; k<=n; k++) {
  780.                 if(tree1[i][k] == tree2[j][k]) nhits1++;
  781.                 if(tree1[i][k] != tree2[j][k]) nhits2++;
  782.             }
  783.             if((nhits1 == last_seq-first_seq+1) || (nhits2 == last_seq-first_seq+1)) hits[i]++;
  784.         }
  785.     }
  786. }
  787.  
  788.  
  789. void print_phylip_tree(char **tree_description, FILE *tree, Boolean bootstrap)
  790. {
  791.     fprintf(tree,"(\n");
  792.  
  793.     two_way_split(tree_description, tree, last_seq-first_seq+1-2,1,bootstrap);
  794.     fprintf(tree,":%7.5f,\n",left_branch[last_seq-first_seq+1-2]);
  795.     two_way_split(tree_description, tree, last_seq-first_seq+1-2,2,bootstrap);
  796.     fprintf(tree,":%7.5f,\n",left_branch[last_seq-first_seq+1-1]);
  797.     two_way_split(tree_description, tree, last_seq-first_seq+1-2,3,bootstrap);
  798.  
  799.     fprintf(tree,":%7.5f);\n",left_branch[last_seq-first_seq+1]);
  800. }
  801.  
  802.  
  803. void two_way_split
  804. (char **tree_description, FILE *tree, int start_row, int flag, Boolean bootstrap)
  805. {
  806.     int row, row2, new_row, col, col2, test_col;
  807.     Boolean single_seq;
  808.  
  809.     if(start_row != last_seq-first_seq+1-2) fprintf(tree,"(\n"); 
  810.  
  811.     for(col=1; col<=last_seq-first_seq+1; col++) {
  812.         if(tree_description[start_row][col] == flag) {
  813.             test_col = col;
  814.             break;
  815.         }
  816.     }
  817.  
  818.     single_seq = TRUE;
  819.     for(row=start_row-1; row>=1; row--) 
  820.         if(tree_description[row][test_col] == 1) {
  821.             single_seq = FALSE;
  822.             new_row = row;
  823.             break;
  824.         }
  825.  
  826.     if(single_seq) {
  827.         tree_description[start_row][test_col] = 0;
  828.         fprintf(tree,"%.10s",names[test_col+first_seq-1]);
  829.     }
  830.     else {
  831.         for(col=1; col<=last_seq-first_seq+1; col++) {
  832.             if((tree_description[start_row][col]==1)&&
  833.                (tree_description[new_row][col]==1))
  834.                 tree_description[start_row][col] = 0;
  835.         }
  836.         two_way_split(tree_description, tree, new_row, 1, bootstrap);
  837.     }
  838.  
  839.     if(start_row == last_seq-first_seq+1-2) return;
  840.  
  841.     fprintf(tree,":%7.5f,\n",left_branch[start_row]);
  842.  
  843.     for(col=1; col<=last_seq-first_seq+1; col++) 
  844.         if(tree_description[start_row][col] == flag) {
  845.             test_col = col;
  846.             break;
  847.         }
  848.     
  849.     single_seq = TRUE;
  850.     for(row=start_row-1; row>=1; row--) 
  851.         if(tree_description[row][test_col] == 1) {
  852.             single_seq = FALSE;
  853.             new_row = row;
  854.             break;
  855.         }
  856.  
  857.     if(single_seq) {
  858.         tree_description[start_row][test_col] = 0;
  859.         fprintf(tree,"%.10s",names[test_col+first_seq-1]);
  860.     }
  861.     else {
  862.         for(col=1; col<=last_seq-first_seq+1; col++) {
  863.             if((tree_description[start_row][col]==1)&&
  864.                (tree_description[new_row][col]==1))
  865.                 tree_description[start_row][col] = 0;
  866.         }
  867.         two_way_split(tree_description, tree, new_row, 1, bootstrap);
  868.     }
  869.  
  870.     fprintf(tree,":%7.5f)\n",right_branch[start_row]);
  871.     if (bootstrap && (boot_totals[start_row]>0))
  872.         fprintf(tree,"%d",(pint)boot_totals[start_row]);
  873. }
  874.  
  875.  
  876.  
  877. void print_tree(char **tree_description, FILE *tree, int *totals)
  878. {
  879.     int row,col;
  880.  
  881.     fprintf(tree,"\n");
  882.  
  883.     for(row=1; row<=last_seq-first_seq+1-3; row++)  {
  884.         fprintf(tree," \n");
  885.         for(col=1; col<=last_seq-first_seq+1; col++) { 
  886.             if(tree_description[row][col] == 0)
  887.                 fprintf(tree,"*");
  888.             else
  889.                 fprintf(tree,".");
  890.         }
  891.         if(totals[row] > 0)
  892.             fprintf(tree,"%7d",(pint)totals[row]);
  893.     }
  894.     fprintf(tree," \n");
  895.     for(col=1; col<=last_seq-first_seq+1; col++) 
  896.         fprintf(tree,"%1d",(pint)tree_description[last_seq-first_seq+1-2][col]);
  897.     fprintf(tree,"\n");
  898. }
  899.  
  900.  
  901.  
  902. int dna_distance_matrix(FILE *tree)
  903. {   
  904.     int m,n,j,i;
  905.     int res1, res2;
  906.         int overspill = 0;
  907.     double p,q,e,a,b,k;    
  908.  
  909.     tree_gap_delete();  /* flag positions with gaps (tree_gaps[i] = 1 ) */
  910.     
  911.     if(verbose) {
  912.         fprintf(tree,"\n");
  913.         fprintf(tree,"\n DIST   = percentage divergence (/100)");
  914.         fprintf(tree,"\n p      = rate of transition (A <-> G; C <-> T)");
  915.         fprintf(tree,"\n q      = rate of transversion");
  916.         fprintf(tree,"\n Length = number of sites used in comparison");
  917.         fprintf(tree,"\n");
  918.         if(tossgaps) {
  919.         fprintf(tree,"\n All sites with gaps (in any sequence) deleted!");
  920.         fprintf(tree,"\n");
  921.         }
  922.         if(kimura) {
  923.         fprintf(tree,"\n Distances corrected by Kimura's 2 parameter model:");
  924.         fprintf(tree,"\n\n Kimura, M. (1980)");
  925.         fprintf(tree," A simple method for estimating evolutionary ");
  926.         fprintf(tree,"rates of base");
  927.         fprintf(tree,"\n substitutions through comparative studies of ");
  928.         fprintf(tree,"nucleotide sequences.");
  929.         fprintf(tree,"\n J. Mol. Evol., 16, 111-120.");
  930.         fprintf(tree,"\n\n");
  931.         }
  932.     }
  933.  
  934.     for(m=1;   m<last_seq-first_seq+1;  ++m)     /* for every pair of sequence */
  935.     for(n=m+1; n<=last_seq-first_seq+1; ++n) {
  936.         p = q = e = 0.0;
  937.         tmat[m][n] = tmat[n][m] = 0.0;
  938.         for(i=1; i<=seqlen_array[first_seq]; ++i) {
  939.             j = boot_positions[i];
  940.                         if(tossgaps && (tree_gaps[j] > 0) ) 
  941.                 goto skip;          /* gap position */
  942.             res1 = seq_array[m+first_seq-1][j];
  943.             res2 = seq_array[n+first_seq-1][j];
  944.             if( (res1 < 0) || (res1 > 4) ||
  945.                             (res2 < 0) || (res2 > 4))  
  946.                 goto skip;          /* gap in a seq*/
  947.             e = e + 1.0;
  948.                         if(res1 != res2) {
  949.                 if(transition(res1,res2))
  950.                     p = p + 1.0;
  951.                 else
  952.                     q = q + 1.0;
  953.             }
  954.                 skip:;
  955.         }
  956.  
  957.  
  958.     /* Kimura's 2 parameter correction for multiple substitutions */
  959.  
  960.         if(!kimura) {
  961.             if (e == 0) {
  962.                 fprintf(stdout,"\n WARNING: sequences %d and %d are non-overlapping\n",m,n);
  963.                 k = 0.0;
  964.                 p = 0.0;
  965.                 q = 0.0;
  966.             }
  967.             else {
  968.                 k = (p+q)/e;
  969.                 if(p > 0.0)
  970.                     p = p/e;
  971.                 else
  972.                     p = 0.0;
  973.                 if(q > 0.0)
  974.                     q = q/e;
  975.                 else
  976.                     q = 0.0;
  977.             }
  978.             tmat[m][n] = tmat[n][m] = k;
  979.             if(verbose)                    /* if screen output */
  980.                 fprintf(tree,        
  981.           "%4d vs.%4d:  DIST = %7.4f; p = %6.4f; q = %6.4f; length = %6.0f\n"
  982.                              ,(pint)m,(pint)n,k,p,q,e);
  983.         }
  984.         else {
  985.             if (e == 0) {
  986.                 fprintf(stdout,"\n WARNING: sequences %d and %d are non-overlapping\n",m,n);
  987.                 p = 0.0;
  988.                 q = 0.0;
  989.             }
  990.             else {
  991.                 if(p > 0.0)
  992.                     p = p/e;
  993.                 else
  994.                     p = 0.0;
  995.                 if(q > 0.0)
  996.                     q = q/e;
  997.                 else
  998.                     q = 0.0;
  999.             }
  1000.  
  1001.             if( ((2.0*p)+q) == 1.0 )
  1002.                 a = 0.0;
  1003.             else
  1004.                 a = 1.0/(1.0-(2.0*p)-q);
  1005.  
  1006.             if( q == 0.5 )
  1007.                 b = 0.0;
  1008.             else
  1009.                 b = 1.0/(1.0-(2.0*q));
  1010.  
  1011. /* watch for values going off the scale for the correction. */
  1012.             if( (a<=0.0) || (b<=0.0) ) {
  1013.                 overspill++;
  1014.                 k = 3.5;  /* arbitrary high score */ 
  1015.             }
  1016.             else 
  1017.                 k = 0.5*log(a) + 0.25*log(b);
  1018.             tmat[m][n] = tmat[n][m] = k;
  1019.             if(verbose)                      /* if screen output */
  1020.                    fprintf(tree,
  1021.              "%4d vs.%4d:  DIST = %7.4f; p = %6.4f; q = %6.4f; length = %6.0f\n"
  1022.                             ,(pint)m,(pint)n,k,p,q,e);
  1023.  
  1024.         }
  1025.     }
  1026.     return overspill;    /* return the number of off-scale values */
  1027. }
  1028.  
  1029.  
  1030. int prot_distance_matrix(FILE *tree)
  1031. {
  1032.     int m,n,j,i;
  1033.     int res1, res2;
  1034.         int overspill = 0;
  1035.     double p,e,a,b,k, table_entry;    
  1036.  
  1037.  
  1038.     tree_gap_delete();  /* flag positions with gaps (tree_gaps[i] = 1 ) */
  1039.     
  1040.     if(verbose) {
  1041.         fprintf(tree,"\n");
  1042.         fprintf(tree,"\n DIST   = percentage divergence (/100)");
  1043.         fprintf(tree,"\n Length = number of sites used in comparison");
  1044.         fprintf(tree,"\n\n");
  1045.             if(tossgaps) {
  1046.             fprintf(tree,"\n All sites with gaps (in any sequence) deleted");
  1047.             fprintf(tree,"\n");
  1048.         }
  1049.             if(kimura) {
  1050.             fprintf(tree,"\n Distances up tp 0.75 corrected by Kimura's empirical method:");
  1051.             fprintf(tree,"\n\n Kimura, M. (1983)");
  1052.             fprintf(tree," The Neutral Theory of Molecular Evolution.");
  1053.             fprintf(tree,"\n Page 75. Cambridge University Press, Cambridge, England.");
  1054.             fprintf(tree,"\n\n");
  1055.             }
  1056.     }
  1057.  
  1058.     for(m=1;   m<nseqs;  ++m)     /* for every pair of sequence */
  1059.     for(n=m+1; n<=nseqs; ++n) {
  1060.         p = e = 0.0;
  1061.         tmat[m][n] = tmat[n][m] = 0.0;
  1062.         for(i=1; i<=seqlen_array[1]; ++i) {
  1063.             j = boot_positions[i];
  1064.                     if(tossgaps && (tree_gaps[j] > 0) ) goto skip; /* gap position */
  1065.             res1 = seq_array[m][j];
  1066.             res2 = seq_array[n][j];
  1067.             if( (res1 < 0)     || (res2 < 0) ||
  1068.                             (res1 > max_aa) || (res2 > max_aa)) 
  1069.                                     goto skip;   /* gap in a seq*/
  1070.             e = e + 1.0;
  1071.                         if(res1 != res2) p = p + 1.0;
  1072.                 skip:;
  1073.         }
  1074.  
  1075.         if(p <= 0.0) 
  1076.             k = 0.0;
  1077.         else
  1078.             k = p/e;
  1079.  
  1080. /* DES debug */
  1081. /* fprintf(stdout,"Seq1=%4d Seq2=%4d  k =%7.4f \n",(pint)m,(pint)n,k); */
  1082. /* DES debug */
  1083.  
  1084.         if(kimura) {
  1085.             if(k < 0.75) { /* use Kimura's formula */
  1086.                 if(k > 0.0) k = - log(1.0 - k - (k * k/5.0) );
  1087.             }
  1088.             else {
  1089.                 if(k > 0.930) {
  1090.                    overspill++;
  1091.                    k = 10.0; /* arbitrarily set to 1000% */
  1092.                 }
  1093.                 else {
  1094.                    table_entry = (k*1000.0) - 750.0;
  1095.                                    k = (double)dayhoff_pams[(int)table_entry];
  1096.                                    k = k/100.0;
  1097.                 }
  1098.             }
  1099.         }
  1100.  
  1101.         tmat[m][n] = tmat[n][m] = k;
  1102.             if(verbose)                    /* if screen output */
  1103.             fprintf(tree,        
  1104.                       "%4d vs.%4d  DIST = %6.4f;  length = %6.0f\n",
  1105.                       (pint)m,(pint)n,k,e);
  1106.     }
  1107.     return overspill;
  1108. }
  1109.  
  1110.  
  1111. void guide_tree(void)
  1112. /* 
  1113.    Routine for producing unrooted NJ trees from seperately aligned
  1114.    pairwise distances.  This produces the GUIDE DENDROGRAMS in
  1115.    PHYLIP format.
  1116. */
  1117. {       char path[FILENAMELEN+1];
  1118.         int i, j;
  1119.         static char **standard_tree;
  1120.  
  1121.         verbose = FALSE;
  1122.   
  1123.         if(empty) {
  1124.                 error("You must load an alignment first");
  1125.                 return;
  1126.         }
  1127.  
  1128.         get_path(seqname,path);
  1129.  
  1130.         standard_tree   = (char **) ckalloc( (nseqs+1) * sizeof (char *) );
  1131.         for(i=0; i<nseqs+1; i++)
  1132.                 standard_tree[i]  = (char *) ckalloc( (nseqs+1) * sizeof(char));
  1133.  
  1134.         nj_tree(standard_tree,clustal_phy_tree_file);
  1135.  
  1136.         print_phylip_tree(standard_tree,phylip_phy_tree_file,FALSE);
  1137.  
  1138.         ckfree((void *)left_branch);
  1139.         ckfree((void *)right_branch);
  1140.         ckfree((void *)kill);
  1141.         ckfree((void *)av);
  1142.         for (i=1;i<nseqs+1;i++)
  1143.                 ckfree((void *)standard_tree[i]);
  1144.         ckfree((void *)standard_tree);
  1145.  
  1146.         fclose(phylip_phy_tree_file);
  1147.  
  1148.         fprintf(stdout,"\nGuide tree        file created:   [%s]\n",
  1149.         phylip_phy_tree_name);
  1150. }
  1151.  
  1152.